Plotting multispectral data¶

Multispectral data can be plotted as:

  1. Individual bands
  2. Spectral indices
  3. True color 3-band images
  4. False color 3-band images

Spectral indices and false color images can both be used to enhance images to clearly show things that might be hidden from a true color image, such as vegetation health.

In [1]:
%store -r band_dict ndvi_da den_redlining_gdf
Try It: Import libraries

Add missing libraries to the imports

In [2]:
import cartopy.crs as ccrs # CRSs
# Interactive tabular and vector data
import hvplot.xarray # Interactive raster
# Overlay plots
import numpy as np # Adjust images
import xarray as xr # Adjust images

import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
import hvplot.pandas

There are many different ways to represent geospatial coordinates, either spherically or on a flat map. These different systems are called Coordinate Reference Systems.

Try It: Prepare to plot

To make interactive geospatial plots, at the moment we need everything to be in the Mercator CRS.

  1. Reproject your area of interest with .to_crs(ccrs.Mercator())
  2. Reproject your NDVI and band raster data using `.rio.reproject(ccrs.Mercator())
Try It: Plot raster with overlay with xarray

Plotting raster and vector data together using pandas and xarray requires the matplotlib.pyplot library to access some plot layour tools. Using the code below as a starting point, you can play around with adding:

  1. Labels and titles
  2. Different colors with cmap and edgecolor
  3. Different line thickness with line_width

See if you can also figure out what vmin, robust, and the .set() methods do.

In [3]:
import xarray as xr

#repjecting
redlining_plot_gdf = den_redlining_gdf.to_crs(ccrs.Mercator())

#reproject raster data with raster in and out
ndvi_plot_da = ndvi_da.rio.reproject(ccrs.Mercator())

#reproj band dict
band_plot_dict = {
    name: da.rio.reproject(ccrs.Mercator())
    for name, da in band_dict.items()
}
In [4]:
ndvi_plot_da.plot(vmin=0, robust=True)
redlining_plot_gdf.plot(ax=plt.gca(), color='none')
plt.gca().set(
    xlabel='', ylabel='', xticks=[], yticks=[]
)
plt.title("NDVI showing vegetation health in Denver", fontsize=16, color='green', loc='center')

#no axis values
plt.xlabel("none", fontsize=12, color="blue", labelpad=10)  #  x-axis label
plt.ylabel("none", fontsize=12, color="green", labelpad=10)  #  y-axis label

plt.legend()

#desc
# plt.text(2, 20, "NDVI values dark is healthy vegetation and green to yellow is unhealthy to no vegetation.", fontsize=10, color="blue")

plt.show()

/tmp/ipykernel_1626/2217805051.py:12: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
  plt.legend()
/opt/conda/lib/python3.11/site-packages/IPython/core/pylabtools.py:170: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  fig.canvas.print_figure(bytes_io, **kw)
No description has been provided for this image

vegetation along river shows 0.0 - 0.1 healthy vegeation, indication of healthy vegetation along stream. The yellow and green is within city indicating little to no vegetation within the city of denver due to building and other structures

Try It: Plot raster with overlay with hvplot

Now, do the same with hvplot. Note that some parameter names are the same and some are different. Do you notice any physical lines in the NDVI data that line up with the redlining boundaries?

In [5]:
(
    ndvi_plot_da.hvplot(
        geo=True,
        xaxis=None, yaxis=None
    )
    * redlining_plot_gdf.hvplot(
        geo=True, crs=ccrs.Mercator(),
        fill_color=None)
)

/opt/conda/lib/python3.11/site-packages/dask/dataframe/__init__.py:42: FutureWarning: 
Dask dataframe query planning is disabled because dask-expr is not installed.

You can install it with `pip install dask[dataframe]` or `conda install dask`.
This will raise in a future version.

  warnings.warn(msg, FutureWarning)
Out[5]:
Try It: Plot bands with linked subplots

The following code will make a three panel plot with Red, NIR, and Green bands. Why do you think we aren’t using the green band to look at vegetation?

In [6]:
raster_kwargs = dict(
    geo=True, robust=True, 
    xaxis=None, yaxis=None
)
(
    (
        band_plot_dict['red'].hvplot(
            cmap='Reds', title='Red Reflectance', **raster_kwargs)
        + band_plot_dict['nir'].hvplot(
            cmap='Greys', title='NIR Reflectance', **raster_kwargs)
        + band_plot_dict['green'].hvplot(
            cmap='Greens', title='Green Reflectance', **raster_kwargs)
    )
    * redlining_plot_gdf.hvplot(
        geo=True, crs=ccrs.Mercator(),
        fill_color=None)
)
Out[6]:

Red is showing red reflectance which represents unhealthy vegetation since healthy vegetation absorbs redlight on the spectrum, and thus will not reflect red. The red indicates unhealthy to no vegeation. The NIR Reflectance should show healthy vegeation where NIR is high reflectance. Along the river should have high NIR reflectance, however, shows low reflectance and not sure why. This also is the case with the green and red refelectance. Along the river should show low red reflectance but in this image shows high red reflectance. Green indicates healthy vegetation since it is being reflected back in thisimage. You can see the healthy vegetation in along the river.

Try It: Plot RBG

The following code will plot an RGB image using both matplotlib and hvplot. It also performs an action called “Contrast stretching”, and brightens the image.

  1. Read through the stretch_rgb function, and fill out the docstring with the rest of the parameters and your own descriptions. You can ask ChatGPT or another LLM to help you read the code if needed! Please use the numpy style of docstrings
  2. Adjust the low, high, and brighten numbers until you are satisfied with the image. You can also ask ChatGPT to help you figure out what adjustments to make by describing or uploading an image.
In [ ]:
rgb_da = (
    xr.concat(
        [
            band_plot_dict['red'],
            band_plot_dict['green'],
            band_plot_dict['blue']
        ],
        dim='rgb')
)

def stretch_rgb(rgb_da, low, high, brighten):
    """
    Short description: 
    True color Satelite image of the city of Denver. 
    Long description...

    Parameters
    ----------
    rgb_da: array-like
      Where the values are set to show contrast low (1), high(99) and brightest(.2)
    param2: ...
      ...

    Returns
    -------
    rgb_da: array-like
      return rgb_da sends the value to produce the output satelite real color image below. 
    """
    p_low, p_high = np.nanpercentile(rgb_da, (low, high))
    rgb_da = (rgb_da - p_low)  / (p_high - p_low) + brighten
    rgb_da = rgb_da.clip(0, 1)
    return rgb_da

rgb_da = stretch_rgb(rgb_da, 1, 99, .2)

rgb_da.plot.imshow(rgb='rgb')
rgb_da.hvplot.rgb(geo=True, x='x', y='y', bands='rgb')
Out[ ]:
No description has been provided for this image
In [8]:
help(stretch_rgb)

Help on function stretch_rgb in module __main__:

stretch_rgb(rgb_da, low, high, brighten)
    Short description: 
    True color Satelite image of the city of Denver. 
    Long description...
    
    Parameters
    ----------
    rgb_da: array-like
      Where the values are set to show contrast low (1), high(99) and brightest(.2)
    param2: ...
      ...
    
    Returns
    -------
    rgb_da: array-like
      return rgb_da sends the value to produce the output satelite real color image below.
Try It: Plot CIR

Now, plot a false color RGB image. CIR images have the following bands:

  • red becomes NIR
  • green becomes red
  • blue becomes green
Looking for an Extra Challenge?: Adjust the levels

You may notice that the NIR band in this image is very bright. Can you adjust it so it is balanced more effectively by the other bands?

In [9]:
rgb_da = (
    xr.concat(
        [
            band_plot_dict['nir'],
            band_plot_dict['red'],
            band_plot_dict['green']
        ],
        dim='rgb')
)

def stretch_rgb(rgb_da, low, high, brighten):
    """
    Short description: 
    false color Satelite image of the city of Denver. 
    Long description...

    Parameters
    ----------
    rgb_da: array-like
      Where the values are set to show contrast low (1), high(99) and brightest(.2)
    param2: ...
      ...

    Returns
    -------
    rgb_da: array-like
      return rgb_da sends the value to produce the output satelite false color image below. 
    """
    p_low, p_high = np.nanpercentile(rgb_da, (low, high))
    rgb_da = (rgb_da - p_low)  / (p_high - p_low) + brighten
    rgb_da = rgb_da.clip(0, 1)
    return rgb_da

rgb_da = stretch_rgb(rgb_da, 1, 99, .2)

rgb_da.plot.imshow(rgb='rgb')
rgb_da.hvplot.rgb(geo=True, x='x', y='y', bands='rgb')
Out[9]:
No description has been provided for this image
In [10]:
help(stretch_rgb)

Help on function stretch_rgb in module __main__:

stretch_rgb(rgb_da, low, high, brighten)
    Short description: 
    false color Satelite image of the city of Denver. 
    Long description...
    
    Parameters
    ----------
    rgb_da: array-like
      Where the values are set to show contrast low (1), high(99) and brightest(.2)
    param2: ...
      ...
    
    Returns
    -------
    rgb_da: array-like
      return rgb_da sends the value to produce the output satelite false color image below.
In [ ]:
%store den_redlining_gdf ndvi_da

Stored 'den_redlining_gdf' (GeoDataFrame)
Stored 'ndvi_da' (DataArray)
In [12]:
%%capture
%%bash
jupyter nbconvert redlining-42-tree-model.ipynb --to html

---------------------------------------------------------------------------
CalledProcessError                        Traceback (most recent call last)
Cell In[12], line 1
----> 1 get_ipython().run_cell_magic('bash', '', 'jupyter nbconvert redlining-34-spectral-plot.ipynb --to html\n')

File /opt/conda/lib/python3.11/site-packages/IPython/core/interactiveshell.py:2541, in InteractiveShell.run_cell_magic(self, magic_name, line, cell)
   2539 with self.builtin_trap:
   2540     args = (magic_arg_s, cell)
-> 2541     result = fn(*args, **kwargs)
   2543 # The code below prevents the output from being displayed
   2544 # when using magics with decorator @output_can_be_silenced
   2545 # when the last Python token in the expression is a ';'.
   2546 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False):

File /opt/conda/lib/python3.11/site-packages/IPython/core/magics/script.py:155, in ScriptMagics._make_script_magic.<locals>.named_script_magic(line, cell)
    153 else:
    154     line = script
--> 155 return self.shebang(line, cell)

File /opt/conda/lib/python3.11/site-packages/IPython/core/magics/script.py:315, in ScriptMagics.shebang(self, line, cell)
    310 if args.raise_error and p.returncode != 0:
    311     # If we get here and p.returncode is still None, we must have
    312     # killed it but not yet seen its return code. We don't wait for it,
    313     # in case it's stuck in uninterruptible sleep. -9 = SIGKILL
    314     rc = p.returncode or -9
--> 315     raise CalledProcessError(rc, cell)

CalledProcessError: Command 'b'jupyter nbconvert redlining-34-spectral-plot.ipynb --to html\n'' returned non-zero exit status 255.